#include "Multigrid.h"

int main()
{
    ///一维possion方程
    int pow_n = 10;
    size_t DIM1_n = 1 << pow_n;
    int dim1_n_layer = pow_n -1;
    auto DIM1_u = [](double x) -> double { return -x*x/2;};
    auto DIM1_f = [](double x) -> double { return 1;};
    double DIM1_x_l = 0.0;
    double DIM1_x_r = 1.0;
    Interator_Methon _type = SOR ;
    double w1 = 2.0 / 3;
    DIM1_Multigrid solution1(dim1_n_layer);
    solution1.set_interator_type(_type, w1);
    solution1.Multigrid_Solution_Dim1(DIM1_x_l, DIM1_x_r, pow_n, dim1_n_layer, DIM1_f, DIM1_u);
    ////直接迭代求解

    
    double h = (DIM1_x_r - DIM1_x_l) / DIM1_n ;
    FDM_DIM1 Dim1_fine(DIM1_n, DIM1_x_l, DIM1_x_r, DIM1_f, DIM1_u);
    Dim1_fine.build();
    SpMat Dim1_fine_mat = Dim1_fine.get_matrix();
    VectorXd Dim1_fine_rhs = Dim1_fine.get_Frhs();
    set<size_t> Dim1_fine_bound = Dim1_fine.get_boundary();
    DealWithBoundary(Dim1_fine_mat, Dim1_fine_rhs, Dim1_fine_bound);
    size_t it1_num = 0;
    VectorXd u_dim1 = VectorXd::Zero(Dim1_fine_rhs.size());
    VectorXd udim1_init = u_dim1;
    VectorXd udim1_realval = VectorXd::Zero(Dim1_fine_rhs.size());
    for(size_t i =0; i <= DIM1_n; i++)
    {
        udim1_realval[i] = DIM1_u(DIM1_x_l + i*h);
    }
    auto t_start = std::chrono::high_resolution_clock::now();
    Interator(Dim1_fine_mat, u_dim1, Dim1_fine_rhs, _type, w1) ;
    it1_num ++; 
    while(it1_num < 1000000 && L2_Error(u_dim1,udim1_init) > 1e-8){
        udim1_init = u_dim1;
        Interator(Dim1_fine_mat, u_dim1, Dim1_fine_rhs, _type, w1) ;
        it1_num ++;    
    }
    auto t_end = std::chrono::high_resolution_clock::now();
    float t_total = std::chrono::duration<float, std::milli>(t_end - t_start).count();
    cout <<"\n\nWithout Multigrid: \n";
    std::cout << "It takes " << t_total <<".ms" << endl;
    cout << "For Dimesion 1:\n" << "The number of interation without multigrid is \n";
    cout << it1_num << endl << "The L2_error of f and Au is :";
    cout << L2_Error( Dim1_fine_rhs ,Dim1_fine_mat*u_dim1) << endl;
    cout << "The L2_error of solution and real_solution is : " ;
    cout << L2_Error(u_dim1, udim1_realval) << endl;
    
    
    
    ///二维possion方程
    /* vector<double> DIM2_lb = { -1, -1};
    vector<double> DIM2_rt = {1, 1};
    int nx_pow = 6;
    int ny_pow = 6;
    int dim2_n_layer = nx_pow -1;
    int acc = 4 ;
    Interator_Methon _type2 = SSOR ;
    double w2 = 2.0/3;
    auto DIM2_u = [](double x, double y) -> double{ return -(x*x + y*y)/4;};
    auto DIM2_f = [](double x, double y) -> double{ return 1; };
    DIM2_Multigrid solution2(dim2_n_layer);
    solution2.set_interator_type(_type2, w2);
    solution2.Multigrid_Solution_Dim2(DIM2_lb, DIM2_rt, nx_pow, ny_pow, dim2_n_layer, DIM2_f, DIM2_u, acc);
    ////直接迭代求解

    size_t DIM2_nx = 1<<nx_pow ;
    size_t DIM2_ny = 1<<ny_pow ;
    
    
    Q1_Rectangel_fem_solution Dim2_fine(DIM2_lb, DIM2_rt, DIM2_nx, DIM2_ny, acc, DIM2_f, DIM2_u);
    Dim2_fine.build();
    SpMat Dim2_fine_mat = Dim2_fine.get_Mat();
    VectorXd Dim2_fine_rhs = Dim2_fine.get_Rhs();
    set<size_t> Dim2_fine_boundary = Dim2_fine.get_boundary();
    DealWithBoundary(Dim2_fine_mat, Dim2_fine_rhs, Dim2_fine_boundary);
    size_t it2_num = 0;
    VectorXd u_dim2 = VectorXd::Zero(Dim2_fine_rhs.size());
    VectorXd u_dim2_init = u_dim2;
    VectorXd u_dim2_realval = VectorXd::Zero(Dim2_fine_rhs.size());
    double DIM2_hx = (DIM2_rt[0] - DIM2_lb[0]) / DIM2_nx;
    double DIM2_hy = (DIM2_rt[1] - DIM2_lb[1]) / DIM2_ny;
    for(size_t i =0; i <= DIM2_ny; i++)
        for(size_t j =0; j <= DIM2_nx; j++)
        {
            u_dim2_realval[i*(DIM2_nx + 1) + j] = DIM2_u(i*DIM2_hy +DIM2_lb[0], j*DIM2_hx +DIM2_lb[1]);
        }
    auto t_start2 = std::chrono::high_resolution_clock::now();
    Interator(Dim2_fine_mat, u_dim2, Dim2_fine_rhs, _type2, w2) ;
    it2_num ++;
    while(it2_num < 1000000 && L2_Error(u_dim2,u_dim2_init) > 1e-8){
        if(it2_num % 10000 == 0)
            cout <<L2_Error(u_dim2,u_dim2_init) << endl;
        u_dim2_init = u_dim2;
        Interator(Dim2_fine_mat, u_dim2, Dim2_fine_rhs, _type2, w2) ;
        it2_num ++;
    }
    auto t_end2 = std::chrono::high_resolution_clock::now();
    float t_total2 = std::chrono::duration<float, std::milli>(t_end2 - t_start2).count();
    cout <<"\n\nWithout Multigrid: \n";
    cout << "For Dimesion 2:\n" << "The number of interation without multigrid is \n";
    cout << it2_num << endl <<"The L2_error of f and Au is: ";
    cout << L2_Error( Dim2_fine_rhs ,Dim2_fine_mat*u_dim2) << endl;
    std::cout << "It takes " << t_total2 <<".ms" << endl;
 */
    return 0 ;
}